Study of unconventional superfluid phases and the phase dynamics in spin-orbit 
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We study the phase distribution and its dynamics in spin-orbit coupled two component ultracold 
Bosons for finite size system. Using an inhomogeneous meanfield analysis we demonstrate how phase 
distribution evolves as we tune the spin-orbit coupling 7 and t, the spin-independent hopping. For 
t » 7 we find the homogeneous superfluid phase. As we increase 7, differences in the phases of 
the order parameter grows leading to twisted superfluid phase. For t ~ 7 competing orderings in 
the phase distribution is seen. At large 7 limit a Ferro-Magnetic stripe ordering appears along the 
diagonal. We explain that this is due to the frustration bought in by the spin-orbit interaction. 
Isolated vortex formation is also shown to appear. We also investigate the possible collective modes. 
In deep superfluid regime we derive the equation of motion for the phases following a semi-classical 
approximation. Imaginary frequencies indicating the damped modes are seen to appear and the 
dynamics of lowest normal modes are discussed. 

PACS numbers: 03.75.Lm, 05.30.Jp, 05.30.Rt 
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I. INTRODUCTION 

The recent advancement in optical lattice experiments 
to investigate the idealized strongly correlated many 
body system has initiated a great interest among the 
condensed matter community 1 . Starting from mimicking 
simple tight binding Hamiltonian in a periodic lattice, it 
can now create more complex situations seen in real ma- 
terials. Creation of artificial abelian or non-abelian gauge 
fields, density-density interaction are some of them to 
mention 2,3 . Experimental realization of Mott-Insulator 
to Superfluid transition for ultracold bosons^ in such 
system became a paradigm of itself. Recently there 
has been experimental realisation to simulate tunable 
spin-orbit coupling in neutral bosons in optical lattice 6,7 
. This has been remarkable because it is known that 
for real material spin-orbit coupling is in essential an 
intrinsic 8 properties of the material and could not be 
controlled. The spin-orbit interaction can change the 
physical properties of the system dramatically. In opti- 
cal lattice the spin-orbit coupling is achieved by Raman 
laser induced transitions between the two internal states 
of a neutral bosonic atom. The resulting spin-orbit inter- 
action could be purely Rashbha 9 type or Dresselhaus 10 
type or suitable combination of both. 

The result of such spin orbit interaction has been stud- 
ied extensively recently-ii~~. In the Mott regime it is 
shown to support exotic magnetic textures, such as vor- 
tex crystals and skyrmion latticed— . The signature 
of the Mott-Insulator to Superfluid transition has been 
shown to be associated with precursor peaks in momen- 
tum distributions^ 7 . - — . Various other equilibrium and 
non-equilibrium dynamics has also been analyzed which 
could have interesting experimental signatures 2 ^. Boson 
fractionalisation has also been proposed and formation of 
twisted superfluid phases has been noticed as a result of 
spin-orbit interactio n 20 ' 22 . It may be mentioned that for 
the fermionic case interesting many body dynamics has 



also been observed 23 . 

The Mott-Insulator to Superfluid transition is well cap- 
tured by Bose-Hubbard model 2 ^— . There are already a 
large number of work done to investigate the low energy 
properties of such Bose-Hubbard model2!~— . However 
much of these work was mainly aimed at investigating the 
systems which are thermodynamically large and in weak 
couple regime. In this work we look into the effect of 
spin-orbit interaction of two component bosons in strong 
coupling limit for different finite size systems. We are mo- 
tivated to look into microscopic manifestation of the spin- 
orbit interaction and various ramifications of superfluid 
order parameter for different system size and different pa- 
rameter regime. For this we employ Gutzwiller projected 
inhomogeneous meanfield treatment 3 - - which seems to be 
pertinent for such small system size. We work in the 
strong coupling limit where the Hubbard interaction is 
the highest energy scale of the problem. This limit en- 
ables us to take the number of states in the Gutzwiller 
projected state to be necessary minimal. Below we de- 
scribe our plan of work. 

In section I, we begin by giving a detail analysis of the 
meanfield procedure and obtain the phase diagram for 
MI-SF transition. Following this, we look into the phases 
and magnitude of the SF order parameter in superfluid 
regime. We show that the phases and the magnitude of 
the SF order parameter respond non-trivially as the pa- 
rameters are varied. We find that when t >> 7, the SF 
phase is described by a homogeneous superfluid where 
the magnitude and the phases of the up spins are spa- 
tially uniform. For intermediate values of t and 7 we find 
that the phases and the amplitudes of both the spins are 
inhomogeneous and shows interesting nontrivial pattern. 
Depending on the relative strength it could be superposi- 
tion of local homogeneous phases and patches where the 
phases form a spiralling pattern. For the limit 7 >> t, 
the phases of the order parameter develops a Ferromag- 
netic order along the diagonal direction followed by pe- 
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FIG. 1: In the left panel we have shown the MI-SF transition 
for A = 0.6, Q, = 0.01 in /j, — t plane. The blue line denotes 
the transition for 7 = and the red line denotes 7 = 0.04. 
In the right panel we have plotted the energy density per 
plaquette with various system sizes for different parameter 
values. The various color represents various set of (7, i) which 
are explained in the text. 



riodic modulations of magnitude of the SF order param- 
eter. We explain that this is due to inherent frustration 
brought in by the spin-orbit interaction. 

In section II, we study the fluctuations around the 
meanfield configuration and investigate into the lowest 
possible excitations. In section III, we study the dynam- 
ics of phases inside the deep SF regime. Assuming that 
the phases of the order parameters are the low energy de- 
grees of freedom in this regime, we deduce the Lagrangc- 
Equation of motion for it. We find the normal modes. 
It appears that due to the constrained collective motion 
imaginary frequency appears signifying damped vibra- 
tion. We also look at the nature of lowest normal modes 
of the vibrations. 



II. MEANFIELD STUDY 

As already mentioned in this work we study a spin- 
orbit coupled two component bosons in square lattice. 
The Hamiltonian for such a system can be written as 
H = Hq + H 1 , where Hq and Hi are given byL2°., 



tfo = £ 



■/mi 



Un ia {n ia - 1) + XU} j n i in i 2 



J2 tatiaha, #1 = *7 £ *i*'( ff X 4)*j 



(ij)o 



E (wUv*i - n*<V**< 



(i) 



Here "Fj = (6ji,&i2)- In the above Hamiltonian /i rep- 
resents the chemical potential, J7 is the Zeeman shift 
between the two species, U is the intraspecies interac- 
tion and A is the on site interspecies interaction. For 
the meanfield analysis we take the Gutzwiller variational 
wave function |>F) = Yii iV'i)) where \ipi) is the wave 



function at a given site H\ \ipi) is given by \ipi) = 
J2 m n fm,n\'mn). As we work in a strong coupling limit 
where U is much larger that t and 7, it is sufficient to 
take states upto 2 particle at a given site. The meanfield 
order parameter is defined as, Aj a — (ipi\bi a \ipi) ■ The 
expression for Aj a in terms of f m n,i are given below, 

Aa = /io,i/oo,t + Am/om + v2/20,i/io,i 

A 4 2 = foi,ifoo,i + fll,iflO,i + ^2/02,1/01,1 (2) 

The first part of the Hamiltonian in Eq.([T]) contains 
the on site interactions and we call it H at which given 
by, 



(H a t)i — 



-Mi(|/km| 2 + I/imI 2 + 2|/ 2CM | 2 ) + A/7|/i M | 2 



M 2 (l/oi,i| 2 + I/imI 2 + 2|/ 2,,| 2 ) + 2£/(|/ 2 cm| 2 + |/ 2,,| 2 ) (3) 

A generic term in Hi can be written as b\ bj t p. The 
meanfield decomposition of it is given by, 



*\j>i,fi = A* 



fijP 



^.iabl - A* a A, 



J i0 



30 



(4) 



After we substitute Eq.Q, Eq.gJ in Eq.© and use 
Eq.([2]) we can write the meanfield decomposed Hamil- 
tonian as, 



H = J2xlF l {fi,X,A J . a ,t, 1 ) Xl 



(5) 



where Xi = {foo,i,fio,i,fai,i,fii,i,ha,i,fo2,i)- The prob- 
lem then reduces to diagonalizing the matrix Fi at every 
site self consistently. The Hamiltonian in Eq. ([S] is still 
a coupled problem. We notice that in the presence of 
spin-orbit coupling Aj can not be taken uniform at each 
site for then the spin-orbit interaction contribute noth- 
ing to the total energy. To find the meanfield solution, 
we start from a given random initial distribution of A,; 
at each site i diagonalize the i^(/i, A, Aj . Q , £,7) at each 
site. We then calculate the new set of Aj corresponding 
to the minimum eigenvalue of Fi . The resulting A j 's are 
fed back into Eq.© until Aj becomes equals to Aj at 
each site i. We do this procedure for approximately 10 4 
random configurations and take the configurations of Aj 
which corresponds to the global minima. In the Fig.(UJ) 
left panel, we show the phase diagram for the MI-SF tran- 
sitions. In the right panel of Fig. ([T]), we have plotted the 
energy density per site with the system size for various 
set of parameter . We find that finite size minimiza- 
tion brings significant variations in the energy density 
with the system size. The various color represents var- 
ious set of parameters (7,*). Red represents (0.1,0.02), 
blue represents (0.02,0.04), green represents (0.03,0.04), 
black is for (0.04,0.04), gray is for (0.06,0.04), orange de- 
notes (0.08,0.04) magenta denotes (0.025,0.04) and cyan 
is for (0.035,0.04). This color scheme is maintained for 
all the figures that will be used later. In the following 
we discuss the textures of the order parameter Aj a for 
different values of t and 7. 
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FIG. 2: Magnitude and phase of the order parameter is plot- 
ted at each site. The arrows represent the phases and the color 
represent the magnitude of the order parameter Aj . The up- 
per panel denotes phase and magnitude for Ai and the lower 
panels are for A2. The left panels denotes the result for 7 = 
0.02, t = 0.04. The middle panels are for 7 = 0.025, t = 0.04 
and the right panels are for 7 = 0.03, t = 0.04. 




FIG. 3: The distribution of phases and the order parameter 
as explained in Fig.((2|). Here the left panels represent the 
result for 7 = 0.035, t — 0.04, the middle panels represent 
7 = 0.04, t = 0.04 and the right panels represent 7 = 0.06, t — 
0.04. 



A. Numerical results 

First we discuss the regime when t » 7 followed by 
the regime where t ~ 7. Lastly we discuss the regime 
where 7 >> t. 

1. Meanfield results when t is large compared to 7. 

In Fig.© we present the resulting distributions of 
phases and the magnitude of the order parameter A. 
The arrows represents the phases and the background 
color represents the relative magnitudes of the order pa- 
rameters. The dark color represents greater magnitude. 
The upper panel is for Ai and the lower panel is for 
A 2 . In the extreme left panel the result is shown for 
t = 0.04,7 = 0.02. We find that the distribution of 
phases Ai are ordered and spatially uniform while that 



of A2 is disordered. The magnitudes of Ai shown form 
a two sublattice structure, however there are degenerate 
solutions with spatially uniform magnitude. It is clear 
that the two sublattice structure is the result of spin-orbit 
interaction. Also we have (Ai) >> (A2). The above tex- 
tures is understood easily as for the presence of f2, the 
system is favoring the condensation of species 1 which re- 
sembles the homogeneous superfluid. The middle panel 
of Fig.©, represents the result for t = 0.04, 7 = 0.03. We 
observe that the phases are no longer uniform leading to 
twisted superfluid phase — . We observe the reduction of 
the ordered pattern of Ai and onset of diagonal order- 
ing. The magnitude of Ai are also random. The A2 also 
shows signature of diagonal ordering. The competition 
of ordering along the two diagonal shows the signature of 
large vortices as seen in the moddlelower panal in Fig.©. 



2. t and 7 is comparable 

The phase textures for this regime could be described 
as follows. We find a competition between local ferro 
magnetic alignment for nearest A^'s and the ferromag- 
netic(FM) ordering along the diagonal neighbors. The 
FM ordering for the neighbors results from direct hop- 
ping. Where as the ferromagnetic ordering along the di- 
agonal is due to the spin-orbit coupling as explained in 
next section. In Fig.©, left panel represents the phase 
distribution for 7 = 0.035, t — 0.04, middle panel is for 
7 = 0.04, t = 0.04 and the right panel is for 7 = 0.06. We 
notice that the minimum energy configuration presented 
here is not unique. There are many degenerate configu- 
rations with identical energy. However the quantum fluc- 
tuations would pick the global minima. For examples, in 
Fig.©, we find the onset of density modulations and no 
vertex formations. There are degenerate meanfield solu- 
tions with completely random density distribution with 
isolated vertex formations. 

3. t is small and 7 is large 

In this regime we notice that the phases forms a ferro- 
magnetic alignment along the diagonal. The magnitude 
of the order parameter are also seen to be modulated. In 
Fig.©, we present the distribution for the phases and 
the magnitude of order parameter for 7 = 0.1, t = 0.02. 
We see that ferromagnetic ordering of phases along the 
diagonal is common. While in the left panel FM order- 
ing happens for both the diagonal, for the middle panel it 
happens for only (1,1) direction. In the left panel isolated 
vertex>ii and anti vertex is seen to appear. To understand 
the phase distribution in this regime it may be useful to 
consider an elementary square plaquette and consider the 
meanfield Hamiltonian for it. Let us consider the hop- 
ping of an up spin under spin-orbit coupling via the sites 
i, i + x, i + x + y and i + y in anti-clockwise direction 
as shown in the right upper panel in fig 01 The mean- 
field decomposition put the following constraints on the 
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FIG. 4: The distribution of phases and the order parame- 
ter as explained in Fig. ^ and ©. The left and the middle 
panal is drawn for 7 = 0.1, t — 0.02. These represent the de- 
generate meanfield configuration. We see that the left panal 
contains vertex and antivertex. The green circle contains the 
vertex configurations and the pink circle contains antivertex. 
In the right panel we have shown the spin-orbit coupled hop- 
ping processes for up-spin starting from site i in anti-clockwise 
direction. 



phases, 



i,l ~ "i+x,2 — if) 0i+x,2 — @i+x+y,l 
H+x+y,l ~ "i+y, 2 = 0, #i+t/,2 ~ 0<,1 = 



2' 

7T 

~2 



(6) 



The above set of equations do not have simultaneous so- 
lutions for all the parameters. One may eleminate #i+ x ,2 
(and &i+y.2 ) from the 1st and 2nd (and 3rd and 4th) to 
solve for Oil and 6i+ x + y ^i to obtain that they are equal, 
the numerical outcome seems to conform this. It then 
poses an ill-defined equation for 9i+ Xt2 (and 6i+ v> 2 ) which 
is fixed to minimize the plaquette energy. The ratio of 
average plaquette energy obtained from numerics to that 
obtained by minimizing a single plaquette is 0.94 which 
is satisfactory. In recapitulation we have shown within 
meanfield how the twisted superfluid phase appears as we 
gradually tune the parameter t and 7 for a tight binding 
Hamiltonian. We have shown the onset of density mod- 
ulations and stripe pattern— for the phases as the 7 is 
increased gradually. 



III. FLUCTUATION AROUND THE 
MEANFIELD 

Here we look into the fluctuations around the mean- 
field solutions obtained in the previous section. To take 
into the role of fluctuation we expand Gutzwiller coeffi- 
cients — fmn,i around its saddle point and expand it by 
fmn,i = fmn,i + Sf?nn,i where / mrM represents the equi- 
librium values. After we substitute it in Eq.© we retain 
the terms which are quadratic in 5f mrit i (and its com- 
plex conjugate). The resulting Hamiltonian then could 



H = &H S V 



(7) 



where * = (ipi, ip2, —ipr, --4n) and ip t = (ip ui ,ipdi)- 
Here ip ui = (<j/oo,i #/io,i £/oi,« 8fu,i Sf 2 o,i Sf 02 ,i) and 
ipdi = VCi- K is clear that H$ represents a 127V x 127V 
Hermitian matrix whose eigenvalues and eigenvectors 
represents the collective modes. It may be noted that the 



substitution, Sf„ 



E* 



,*e 



ikr 



-ikr 



v mn ,ke ""' , does 

,,,,,, j s are not trans- 



not simplify the problem as the /, 
lational invariant. We denote the eigenvalue closer to 
absolute zero by E . The E is a measure of possible 
Goldstone modes of the system and is shown in . We 
find that for t » 7, the system always find zero energy 
modes. For t ~ 7, where the phases are disordered we 
also find similar behavior. However for 7 >> t, we find 
that E is ~ 10 4 times larger than the other parameter 
regime. However the Eq scales to lower values monotoni- 
cally as we increase the system size. The gradual decrease 
of Eq with system size TV indicates that it is approach- 
ing to possible zero energy modes. The reason that Eq 
for 7 >> t is larger than other cases by few thousand 
order is the following. For t » 7 the uniform phase 
distribution always find Goldstone modes and there is 
no frustration in the system also. For t ~ 7, the spins 
are disordered and random. Thus it is easily possible to 
re-distribute the phases to have zero energy eigenmodcs 
which is nearly degenerate with the original solutions. 
However for 7 >> t, the distribution of phases and the 
magnitudes are governed by the frustration bought in by 
spin-orbit coupling. The degenerate solutions in this case 
as seen from Fig. ((4]) are not easily connected. Thus the 
collective excitations costs finite energy than the other 
cases. However as we increase the system size, we expect 
that the degenerate solutions are easily obtained from 
one other leading to zero energy mode. We also observe 
that the eigenvalues of the collective modes form three 
distinct bands. This is clear from Eq.©. The fluctua- 
tion of /2,o or /o,2 yields the bands around U. While the 
fluctuation of /n yields the bands around A/2. The fluc- 
tuation of /i : o, /o,i and / ,o constitutes the lower bands. 
We denote these three bands by E 2 , E\ and E respec- 
tively. In the right panel of Fig.©, we have plotted the 
band-width with the system sizes for different parameter 
values. In the left panel of Fig.© we have plotted the 
bandwidth of E\ and the right panel is for Eq. It appears 
that for a given t, the bandwidth is inversely proportional 
to 7. Also more the value of 7, the bandwidth oscillates 
more with the system sizes. We notice that the bands 
E 2 and E\ are symmetric but E is not because of the 
presence of fl. 



IV. DYNAMICS OF THE PHASES 

Now we turn out attention to the deep inside the su- 
perfluid regime where one may neglect the fluctuations 
of the magnitude of the order parameter and consider 





FIG. 5: In the left a measure of zero energy eigenmodes Eq 
due to the collective motions has been shown. In the right 
panel the bandwidth around E2 has been plotted. For both 
the figure horizontal axis represnts the length a N x N lattice. 
The various colors represent various set of parameters as given 
in the text. 
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FIG. 6: In the left panel the bandwidth around E\ has been 
shown. The right panel is for Eo- In both the figure the 
horizontal axis represents the length of & N x N lattice. 



the phases as the only relevant degree of freedom. Fol- 
lowing a semi-classical approximation, we deduce the 
Lagrangian and the equation of motion for the phases 
and determine the normal modes of the vibrations. The 
mcanficld decomposition of Eq.([T]) could be written as, 

H = J~) Ha,i(n a ,i) + —(n a j) 2 +XU(n a ,i)(np >i ) 



(ij) 



(8) 



In the above \ij, a p denotes a general hopping param- 
eter. The main disadvantage of Eq. (jHJ is that all the 
variables commute with each other and bear no signa- 
ture of the original bosonic commutation relations. To 
derive Lagrangian of the phases of the order parameter 
A, we follow the procedure i n 33 ' 34 . Translating the orig- 
inal bosonic commutators to the commutation relations 
of the meanfield variables, we find that, 



[ni, 61] = -61 -> [(m),A 1 ] = -Ai 



(9) 



Writing Ai = e iSl |Ai| and keeping |Ai| constant we ob- 
tain, 



l(ni),e l 



(10) 



Expanding e 101 and keeping only the lowest order term we 
obtain for 8\ — > 0, the following commutation relations, 



[<ni),0i]=*, [(nx) 2 A] 



2i(n 1 ) 



(11) 



The above procedure yields the following coupled equa- 
tions to be solved for the ^f- and 



dt 



dt 
dt 



(n + tt + 


U 
~2 


(m- 


-Q + 


U 
~2 



U(nu) 

U(n 2z ) 



\U{n 2l ) 
XU(n u ) (12) 



Solving for (flu) and (n 2 i) from the above two equations 
and substituting in the Hamiltonian, Eq. [8j we obtain 
the following equations, 



H, 



M( 



d6u r+& 2 



dt 



dt 



Bi 



dt 



B. 



' B3 -df-^ + F{eil ' 9i2) 



B 4 



d6 2l 
dt 

(13) 



Where F(8n,0i 2 ) is given in the appendix. Expressions 
for Aj's are also given in the appendix. To derive the E-L 
equations of motion, we introduce the relative and total 
phase by the relation, 9n = 6 c i+6i r , 9 2 i = Q ci —Q ir After 
inserting the above change of variables we can rewrite Eq 
fTTlas follows, 

i 

+ F(6 lr ,9. lc ) + Y,a icr (14) 

Here Ti/ 2 — 2B ± B 3 Using the above equations, we 
write the resulting Lagrangian and the equation of mo- 
tion below, 

C = 5^Ti(4 + a c ) 2 + T 2 (6„. + a r f - F{0„, 6 lc ) 

i 

■a _ dF(6 ic ,0 ir ) ■■ dF(9 ic ,9 ir ) 



T 1 d6 z , 



T 2 d9 v 



(15) 



In the last equation we have deliberately omitted the in- 
consequential constant term ^ai cr . After simplifying 
the r.h.s of Eqll5] and subsequently expanding upto lin- 
ear term we can rewrite it is, = MQ. Where for a 
system of N x N lattice O is a column matrix with 2N 2 
element such that 6^ = 9i C and 9Ar2 +i = 9 ir where i runs 
from 1 to N 2 . M is a 2 N 2 x 2 A^ 2 matrix. The eigenvalues 
of the matrix M yields the normal modes. We find that 
the due to the presence of 7, the normal modes develop 
negative eigenvalues signifying damped modes. In fig [7] 
we have plotted schematically the lowest normal modes 
for three different regime. In all the plot the blue re- 
gion denotes displacements of phases in forward direction 
(anti-clockwise rotation) and the white regions denotes 
displacements in the backward directions (clockwise ro- 
tation). The right panel denotes the case for 7 >> t, 
the middle panel denotes 7 ~ t and the right panel is for 
t >> 7. In each of these panel the upper one denotes the 
displacement for species 1 and the lower panel describe 
the displacements for species 2. Looking at the upper 
panel we find that for the 7 >> t, there is tendency of 




FIG. 7: We have shown the nature of vibration for the lowest 
normal modes. In each panel the upper panales denotes Q\ 
and the lower panels are for 82- The right panel denotes 
the the case 7 < t, (7 = 0.02, t — 0.04). The middle panal 
reprensents 7 ~ t, (7 = 0.04, t — 0.04). The right panel 
represent 7 > t, (7 = 0.1, t — 0.02). The white region denotes 
the motion in clockwise direction and the blue region denotes 
motion in the anti-clockwise direction. 



phases to move synchronously along the diagonal which 
is expected. However for the middle panel and the left 
panel there is a preferences in horizontal ordering and 



patches of areas vibrating in breathing modes. For the 
species 2, the left panel, we find similar behavior though 
region executing breathing modes are larger. 



V. DISCUSSION 

To summarize we have explored the different phases 
that might occur for a spin-orbit coupled bosons in the 
optical lattice. We have extensively studied the distri- 
bution of phases and the magnitude of order parameter 
for varying finite size system using an inhomogeneous 
meanfield analysis. We have shown that for a given t, 
as we increase the spin-orbit interaction 7, we observe 
the destruction of normal homogeneous superfiuid phase 
and onset of twisted superfiuid phases. At large 7 limit 
an interesting ordering along the diagonal appears. We 
have also investigated the fluctuation around the mean- 
field and shows the existence of Goldstone modes. The 
scaling of minimum energy excitations with system size 
has also been shown. Finally, using semiclassical approx- 
imation we derived the equation of motion for the phases 
and derive the normal modes of vibrations. We think 
that some of the results may have interesting experimen- 
tal signatures in the light of recent experiments. 
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b\ = A (— - a 2 ) + A u(a 2 - — ) +wAA (ai - — ) 
h = A (— - ai) + A u(ai - — -) +MAA (a 2 - — ) 
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u u X 
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2 2 u(A z — 1) 
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